#!/usr/bin/python
# -*- coding: utf-8 -*-

import os
import math

import numpy
import xml.etree.ElementTree as ET

import lib
from lib.nmr import peaks, integrals

class NMR(object):
    s_storagefilename = 'glabsnmr.glne'
    
class NMR_1D(NMR):
    sPath = ''
    
    y_values = []
    
    o_peaks = None
    o_integrals = None
    
    f_integral_reference = None
    f_custom_l_border = 0.0
    f_custom_r_border = 0.0
    
    def __init__(self, path):
        self.sPath = path
        
        self.o_peaks = peaks.PeakList(self)
        self.o_integrals = integrals.IntegralList(self)
        
        self.f_integral_reference = None
        
        self.o_workspace = lib.Workspace._the_instance
    
    def searchPeaks(self, y_axis = None, delta = 0, threshold = None, lookahead = 2, x_axis = None):
        """
        Taken from https://gist.github.com/1178136 with no licence remark. I hope this is okay!
        Modified to fit the special needs of a nmr spectra.
        
        Converted from/based on a MATLAB script at http://billauer.co.il/peakdet.html
        
        Algorithm for detecting local maximas and minmias in a signal.
        Discovers peaks by searching for values which are surrounded by lower
        or larger values for maximas and minimas respectively
        
        keyword arguments:
        y_axis -- A list containg the signal over which to find peaks
        x_axis -- A x-axis whose values correspond to the 'y_axis' list and is used
            in the return to specify the postion of the peaks. If omitted the index
            of the y_axis is used. (default: None)
        lookahead -- (optional) distance to look ahead from a peak candidate to
            determine if it is the actual peak (default: 500) 
            '(sample / period) / f' where '4 >= f >= 1.25' might be a good value
        delta -- (optional) this specifies a minimum difference between a peak and
            the following points, before a peak may be considered a peak. Useful
            to hinder the algorithm from picking up false peaks towards to end of
            the signal. To work well delta should be set to 'delta >= RMSnoise * 5'.
            (default: 0)
                Delta function causes a 20% decrease in speed, when omitted
                Correctly used it can double the speed of the algorithm
                
        threshold -- (optional) this specifies a minimum value which y should have to 
            classify a Peak as a peak. It is normally three times the standard derivation 
            which is considered as the detection limit.
            (default: None)
        
        return -- two lists [maxtab, mintab] containing the positive and negative
            peaks respectively. Each cell of the lists contains a tupple of:
            (position, peak_value) 
            to get the average peak value do 'numpy.mean(maxtab, 0)[1]' on the results
     
        """
        maxtab = []
        mintab = []
        dump = []   #Used to pop the first hit which always if false
        
        if y_axis == None:
            y_axis = self.y_values
        
        length = len(y_axis)
        
        if x_axis is None:
            x_axis = range(length)
        
        #perform some checks
        if length != len(x_axis):
            raise ValueError, "Input vectors y_axis and x_axis must have same length"
        if lookahead < 1:
            raise ValueError, "Lookahead must be above '1' in value"
        if not (numpy.isscalar(delta) and delta >= 0):
            raise ValueError, "delta must be a positive number"
        
        #needs to be a numpy array
        y_axis = numpy.asarray(y_axis)
        
        #maxima and minima candidates are temporarily stored in
        #mx and mn respectively
        mn, mx = numpy.Inf, -numpy.Inf
        
        #Only detect peak if there is 'lookahead' amount of points after it
        for index, (x, y) in enumerate(zip(x_axis[:-lookahead], y_axis[:-lookahead])):
            if threshold != None:
                if y < threshold:
                    continue
            
            if y > mx:
                mx = y
                mxpos = x
            if y < mn:
                mn = y
                mnpos = x
            
            ####look for max####
            if y < mx-delta and mx != numpy.Inf:
                #Maxima peak candidate found
                #look ahead in signal to ensure that this is a peak and not jitter
                if y_axis[index:index+lookahead].max() < mx:
                    maxtab.append((mxpos, mx))
                    dump.append(True)
                    #set algorithm to only find minima now
                    mx = numpy.Inf
                    mn = numpy.Inf
            
            ####look for min####
            if y > mn+delta and mn != -numpy.Inf:
                #Minima peak candidate found 
                #look ahead in signal to ensure that this is a peak and not jitter
                if y_axis[index:index+lookahead].min() > mn:
                    mintab.append((mnpos, mn))
                    dump.append(False)
                    #set algorithm to only find maxima now
                    mn = -numpy.Inf
                    mx = -numpy.Inf
        
        
        #Remove the false hit on the first value of the y_axis
        try:
            if dump[0]:
                maxtab.pop(0)
                #print "pop max"
            else:
                mintab.pop(0)
                #print "pop min"
            del dump
        except IndexError:
            #no peaks were found, should the function return empty lists?
            pass
        
        return maxtab, mintab
        
    def searchPeak(self, i, sharpness = 10, delta = 0, threshold = None):
        sub = int(i-math.ceil(sharpness/2))
        sup = int(i+math.ceil(sharpness/2))
        
        #print "Peak-Lookup in %s..%i..%s" % (sub, i, sup)
        
        x_vals = range(sub, sup, 1)  
        #print self.y_values[sub:sup]      
        maxtab, mintab = self.searchPeaks(self.y_values[sub:sup], delta, threshold, 2, x_axis = x_vals)
        
        max_val = (0, 0)
        i = 0
        j = len(maxtab)
        
        for i in range(0, j):
            if maxtab[i][1] > max_val[1]:
                max_val = maxtab[i]
                
        return max_val        
        
    def searchPeak_old(self, i, sharpness = 10):
        # Schauen, ob schonmal abgeleitet wurde - nicht? Dann ableiten!
        if self.y_values_d1 == {}:
            self.y_values_d1 = NMR_1D.derive(self.y_values)
            
        # Bestimme unteres Ende
        lowerRange = int(i - math.floor(abs(sharpness)/2) - 1)
        if lowerRange < 0:
            lowerRange = 0
        
        # Bestimme oberes Ende
        upperRange = int(i + math.floor(abs(sharpness)/2) + 1)
        if upperRange >= len(y_values):
            upperRange = len(y_values)
        
        # Maximum suchen
        maxX = 0
        maxY = 0
        mindY = 0
        for x in range(lowerRange, upperRange):
            # Muss grösser sein als das bisherige Maximum - und die erste Ableitung kleiner.
            # Sollte so ziemlich alle Peaks finden... hoffentlich.
            if self.y_values[x] > maxY and self.y_values_d1[x] < mindY:
                maxY = self.y_values[y]
                mindY = self.y_values[y]
                maxX = x
        
        # Ignoriere Peaks bei lowerRange und upperRange
        if maxX == lowerRange or maxX == upperRange:
            return False
        else:
            return maxX        

    def addPeak(self, iN):
        if self.o_peaks == None:
            self.o_peaks = peaks.PeakList(self)
        
        self.o_peaks.addPeak(iN)
        
    def getPeaks(self):
        return self.o_peaks
        
    def deletePeak(self, peak = None, n = None):
        self.o_peaks.deletePeak(peak = peak, n = n)
        
    def addIntegral(self, start, end):
        if self.o_integrals == None:
            self.o_integrals = integrals.IntegralList(self)
            
        self.o_integrals.addIntegral(start, end)
        
    def getIntegrals(self):
        return self.o_integrals
        
    def deleteIntegral(self, integral = None, start = None, end = None):
        self.o_integrals.deleteIntegral(integral = integral, start = start, end = end)
        
    def n2height(self, n):
        return self.y_values[n]
        
    def getSlice(self, n_from, n_to):
        return self.y_values[n_from:n_to+1]
    
    def setIntegralReference(self, ref):
        self.f_integral_reference = ref
        
    def getIntegralReference(self):
        return self.f_integral_reference if self.f_integral_reference != None else 100.0
        
    def hasIntegralReference(self):
        return True if self.f_integral_reference != None else False
        
    def save(self):
        """ Save the current nmr! """
        # Root-Element
        nmr = ET.Element("nmr")
        
        printSettings = ET.SubElement(nmr, "print", 
            leftborder = str(self.getCustomLeftBorder()), 
            rightborder = str(self.getCustomRightBorder())
        )
        
        # Create Peak-Elements
        peaks = ET.SubElement(nmr, "peaks")
        if self.o_peaks != None:
            for peak in self.o_peaks:
                p = ET.SubElement(peaks, "peak", n = str(peak.getN()))
            
        # Create Integral-Elements
        if self.hasIntegralReference():
            integrals = ET.SubElement(nmr, "integrals", reference = str(self.f_integral_reference))
        else:
            integrals = ET.SubElement(nmr, "integrals")
            
        if self.o_integrals != None:
            for integral in self.o_integrals:
                i = ET.SubElement(integrals, "integral", start = str(integral.getStart()), end = str(integral.getEnd()))
        
        # Write to disk
        ET.ElementTree(nmr).write(os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'glabsnmrdata'))
        
    def load(self):
        path = os.path.join(self.o_workspace.getNMRDir(), self.sPath, 'glabsnmrdata')
        if os.path.exists(path):
            tree = ET.parse(path)
            nmr = tree.getroot()
            
            # Print Settings
            printSettings = nmr.find("print")
            if printSettings != None:
                leftborder = printSettings.get("leftborder")
                rightborder = printSettings.get("rightborder")
                
                if leftborder != None and rightborder != None:
                    self.setCustomLeftBorder(numpy.float(leftborder))
                    self.setCustomRightBorder(numpy.float(rightborder))
            
            # Peaks
            peaks = nmr.find("peaks")
            for peak in peaks:
                n = peak.get("n")
                if n != None:
                    self.addPeak(int(n))
                    
            # Integrals
            integrals = nmr.find("integrals")
            ref = integrals.get("reference")
            if ref != None:
                self.setIntegralReference(numpy.float(ref))
            
            for integral in integrals:
                start = integral.get("start")
                end = integral.get("end")
                
                if start != None and end != None:
                    self.addIntegral(int(start), int(end))
                
        else:
            return False
